﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Triangulator;

namespace SoilProfiler
{
    public  class Triangulate2
    {
        public static bool   circum(
          double x1, double y1, double x2,
          double y2, double x3, double y3,
          ref double xc, ref double yc, ref double r)
        {
            // Calculation of circumscribed circle coordinates and
            // squared radius

            const double eps = 1e-6;
            const double big = 1e12;
            bool result = true;
            double m1, m2, mx1, mx2, my1, my2, dx, dy;

            if ((Math.Abs(y1 - y2) < eps) && (Math.Abs(y2 - y3) < eps))
            {
                result = false;
                xc = x1; yc = y1; r = big;
            }
            else
            {
                if (Math.Abs(y2 - y1) < eps)
                {
                    m2 = -(x3 - x2) / (y3 - y2);
                    mx2 = (x2 + x3) / 2;
                    my2 = (y2 + y3) / 2;
                    xc = (x2 + x1) / 2;
                    yc = m2 * (xc - mx2) + my2;
                }
                else if (Math.Abs(y3 - y2) < eps)
                {
                    m1 = -(x2 - x1) / (y2 - y1);
                    mx1 = (x1 + x2) / 2;
                    my1 = (y1 + y2) / 2;
                    xc = (x3 + x2) / 2;
                    yc = m1 * (xc - mx1) + my1;
                }
                else
                {
                    m1 = -(x2 - x1) / (y2 - y1);
                    m2 = -(x3 - x2) / (y3 - y2);
                    if (Math.Abs(m1 - m2) < eps)
                    {
                        result = false;
                        xc = x1;
                        yc = y1;
                        r = big;
                    }
                    else
                    {
                        mx1 = (x1 + x2) / 2;
                        mx2 = (x2 + x3) / 2;
                        my1 = (y1 + y2) / 2;
                        my2 = (y2 + y3) / 2;
                        xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
                        yc = m1 * (xc - mx1) + my1;
                    }
                }
            }
            dx = x2 - xc;
            dy = y2 - yc;
            r = dx * dx + dy * dy;
            return result;
        }
    }

    public class DelaunyTriangulation
    {

        /// <param name="Vertex">List of vertices to triangulate.</param>
        /// <returns>Triangles referencing vertex indices arranged in clockwise order</returns>
        public List<Triangulator.Geometry.Triangle> Triangulate(List<Triangulator.Geometry.Point> Vertex)
        {
            int npts = Vertex.Count;
            int nv=npts;
            int trimax = 4 * npts;
            const int maxpoints = 32767;
            // Find the maximum and minimum vertex bounds.
            // This is to allow calculation of the bounding supertriangle
            double xmin = Vertex[0].X;
            double ymin = Vertex[0].Y;
            double xmax = xmin;
            double ymax = ymin;
            for (int ii = 1; ii < nv; ii++)
            {
                if (Vertex[ii].X < xmin) xmin = Vertex[ii].X;
                if (Vertex[ii].X > xmax) xmax = Vertex[ii].X;
                if (Vertex[ii].Y < ymin) ymin = Vertex[ii].Y;
                if (Vertex[ii].Y > ymax) ymax = Vertex[ii].Y;
            }

            double dx = xmax - xmin;
            double dy = ymax - ymin;
            double dmax = (dx > dy) ? dx : dy;

            double xmid = (xmax + xmin) * 0.5;
            double ymid = (ymax + ymin) * 0.5;
            /*
            int i, j, k, ntri, ned, status1 = 0, status2 = 0;
            bool status;
            i = Vertex.Count();
            double[] ptx = new double[maxpoints + 3];
            double[] pty = new double[maxpoints + 3];
            double[] ptz = new double[maxpoints + 3];
            int[] pt1 = new int[maxpoints * 2 + 1];
            int[] pt2 = new int[maxpoints * 2 + 1];
            int[] pt3 = new int[maxpoints * 2 + 1];
            double[] cex = new double[maxpoints * 2 + 1];
            double[] cey = new double[maxpoints * 2 + 1];
            double[] rad = new double[maxpoints * 2 + 1];

            for (int gg = 0; gg < Vertex.Count(); gg++)
			{
			   ptx[gg]=Vertex[gg].X;
               pty[gg]=Vertex[gg].Y;
               ptz[gg]=0;
			}


            int[] ed1 = new int[maxpoints * 2 + 1];
            int[] ed2 = new int[maxpoints * 2 + 1];
            ptx[i] = xmid - (90 * (dx + dy)) - 100;
            pty[i] = ymid - (50 * (dx + dy)) - 100;
            ptz[i] = 0;
            pt1[0] = i;
            i++;
            ptx[i] = xmid + (90 * (dx + dy)) + 100;
            pty[i] = ymid - (50 * (dx + dy)) - 100;
            ptz[i] = 0;
            pt2[0] = i;
            i++;
            ptx[i] = xmid;
            pty[i] = ymid + 100 * (dx + dy + 1);
            ptz[i] = 0;
            pt3[0] = i;
            ntri = 1;
            Triangulate2.circum(
   ptx[pt1[0]], pty[pt1[0]], ptx[pt2[0]],
   pty[pt2[0]], ptx[pt3[0]], pty[pt3[0]],
   ref cex[0], ref cey[0], ref rad[0]
 );


                     for (i = 0; i < Vertex.Count(); i++)
                    {
                      ned = 0;
                      xmin = ptx[i]; ymin = pty[i];
                      j = 0;
                      while (j < ntri)
                      {
                        dx = cex[j] - xmin; dy = cey[j] - ymin;
                        if (((dx * dx) + (dy * dy)) < rad[j])
                        {
                          ed1[ned] = pt1[j]; ed2[ned] = pt2[j];
                          ned++;
                          ed1[ned] = pt2[j]; ed2[ned] = pt3[j];
                          ned++;
                          ed1[ned] = pt3[j]; ed2[ned] = pt1[j];
                          ned++;
                          ntri--;
                          pt1[j] = pt1[ntri];
                          pt2[j] = pt2[ntri];
                          pt3[j] = pt3[ntri];
                          cex[j] = cex[ntri];
                          cey[j] = cey[ntri];
                          rad[j] = rad[ntri];
                          j--;
                        }
                        j++;
                      }

                      for (j = 0; j < ned - 1; j++)
                          for (k = j + 1; k < ned; k++)
                              if ((ed1[j] == ed2[k]) && (ed2[j] == ed1[k]))
                              {
                                  ed1[j] = -1; ed2[j] = -1; ed1[k] = -1; ed2[k] = -1;
                              }

                      for (j = 0; j < ned; j++)
                          if ((ed1[j] >= 0) && (ed2[j] >= 0))
                          {
                              pt1[ntri] = ed1[j]; pt2[ntri] = ed2[j]; pt3[ntri] = i;
                              status =
                                Triangulate2.circum(
                                  ptx[pt1[ntri]], pty[pt1[ntri]], ptx[pt2[ntri]],
                                  pty[pt2[ntri]], ptx[pt3[ntri]], pty[pt3[ntri]],
                                  ref cex[ntri], ref cey[ntri], ref rad[ntri]
                                );
                              if (!status)
                              {
                                  status1++;
                              }
                              ntri++;
                          }
                    }

                     // removal of outer triangles
                     i = 0;
                     while (i < ntri)
                     {
                         if ((pt1[i] >= npts) || (pt2[i] >= npts) || (pt3[i] >= npts))
                         {
                             ntri--;
                             pt1[i] = pt1[ntri];
                             pt2[i] = pt2[ntri];
                             pt3[i] = pt3[ntri];
                             cex[i] = cex[ntri];
                             cey[i] = cey[ntri];
                             rad[i] = rad[ntri];
                             i--;
                         }
                         i++;
                     }

                     List<Triangulator.Geometry.Triangle> tria = new List<Triangulator.Geometry.Triangle>();
                     for (int itt = 0; itt < ntri; itt++)
                     {
                         tria.Add(new Triangulator.Geometry.Triangle(pt1[itt], pt2[itt], pt3[itt]));
                     }

                     return tria;
*/

            // Set up the supertriangle
            // This is a triangle which encompasses all the sample points.
            // The supertriangle coordinates are added to the end of the
            // vertex list. The supertriangle is the first triangle in
            // the triangle list.
            Vertex.Add(new Triangulator.Geometry.Point((xmid - 2 * dmax), (ymid - dmax)));
            Vertex.Add(new Triangulator.Geometry.Point(xmid, (ymid + 2 * dmax)));
            Vertex.Add(new Triangulator.Geometry.Point((xmid + 2 * dmax), (ymid - dmax)));
            List<Triangulator.Geometry.Triangle> Triangle = new List<Triangulator.Geometry.Triangle>();
            Triangle.Add(new Triangulator.Geometry.Triangle(nv, nv + 1, nv + 2)); //SuperTriangle placed at index 0

            // Include each point one at a time into the existing mesh
            for (int i = 0; i < nv; i++)
            {
                List<Triangulator.Geometry.Edge> Edges = new List<Triangulator.Geometry.Edge>(); //[trimax * 3];
                // Set up the edge buffer.
                // If the point (Vertex(i).x,Vertex(i).y) lies inside the circumcircle then the
                // three edges of that triangle are added to the edge buffer and the triangle is removed from list.
                for (int j = 0; j < Triangle.Count; j++)
                {
                    if (InCircle(Vertex[i], Vertex[Triangle[j].p1], Vertex[Triangle[j].p2], Vertex[Triangle[j].p3]))
                    {
                        Edges.Add(new Triangulator.Geometry.Edge(Triangle[j].p1, Triangle[j].p2));
                        Edges.Add(new Triangulator.Geometry.Edge(Triangle[j].p2, Triangle[j].p3));
                        Edges.Add(new Triangulator.Geometry.Edge(Triangle[j].p3, Triangle[j].p1));
                        Triangle.RemoveAt(j);
                        j--;
                    }
                }
                if (i >= nv) continue; //In case we the last duplicate point we removed was the last in the array

                // Remove duplicate edges
                // Note: if all triangles are specified anticlockwise then all
                // interior edges are opposite pointing in direction.
                for (int j = Edges.Count - 2; j >= 0; j--)
                {
                    for (int k = Edges.Count - 1; k >= j + 1; k--)
                    {
                        if (Edges[j].Equals(Edges[k]))
                        {
                            Edges.RemoveAt(k);
                            Edges.RemoveAt(j);
                            k--;
                            continue;
                        }
                    }
                }
                // Form new triangles for the current point
                // Skipping over any tagged edges.
                // All edges are arranged in clockwise order.
                for (int j = 0; j < Edges.Count; j++)
                {
                    if (Triangle.Count >= trimax)
                        throw new ApplicationException("Exceeded maximum edges");
                    Triangle.Add(new Triangulator.Geometry.Triangle(Edges[j].p1, Edges[j].p2, i));
                }
                Edges.Clear();
                Edges = null;
            }
            // Remove triangles with supertriangle vertices
            // These are triangles which have a vertex number greater than nv
            for (int i = Triangle.Count - 1; i >= 0; i--)
            {
                if (Triangle[i].p1 >= nv || Triangle[i].p2 >= nv || Triangle[i].p3 >= nv)
                    Triangle.RemoveAt(i);
            }
            //Remove SuperTriangle vertices
            Vertex.RemoveAt(Vertex.Count - 1);
            Vertex.RemoveAt(Vertex.Count - 1);
            Vertex.RemoveAt(Vertex.Count - 1);
            Triangle.TrimExcess();
            return Triangle;
        }

        /// <summary>
        /// Returns true if the point (p) lies inside the circumcircle made up by points (p1,p2,p3)
        /// </summary>
        /// <remarks>
        /// NOTE: A point on the edge is inside the circumcircle
        /// </remarks>
        /// <param name="p">Point to check</param>
        /// <param name="p1">First point on circle</param>
        /// <param name="p2">Second point on circle</param>
        /// <param name="p3">Third point on circle</param>
        /// <returns>true if p is inside circle</returns>
        private bool InCircle(Triangulator.Geometry.Point p, Triangulator.Geometry.Point p1, Triangulator.Geometry.Point p2, Triangulator.Geometry.Point p3)
        {
            //Return TRUE if the point (xp,yp) lies inside the circumcircle
            //made up by points (x1,y1) (x2,y2) (x3,y3)
            //NOTE: A point on the edge is inside the circumcircle

            if (System.Math.Abs(p1.Y - p2.Y) < double.Epsilon && System.Math.Abs(p2.Y - p3.Y) < double.Epsilon)
            {
                //INCIRCUM - F - Points are coincident !!
                return false;
            }

            double m1, m2;
            double mx1, mx2;
            double my1, my2;
            double xc, yc;

            if (System.Math.Abs(p2.Y - p1.Y) < double.Epsilon)
            {
                m2 = -(p3.X - p2.X) / (p3.Y - p2.Y);
                mx2 = (p2.X + p3.X) * 0.5;
                my2 = (p2.Y + p3.Y) * 0.5;
                //Calculate CircumCircle center (xc,yc)
                xc = (p2.X + p1.X) * 0.5;
                yc = m2 * (xc - mx2) + my2;
            }
            else if (System.Math.Abs(p3.Y - p2.Y) < double.Epsilon)
            {
                m1 = -(p2.X - p1.X) / (p2.Y - p1.Y);
                mx1 = (p1.X + p2.X) * 0.5;
                my1 = (p1.Y + p2.Y) * 0.5;
                //Calculate CircumCircle center (xc,yc)
                xc = (p3.X + p2.X) * 0.5;
                yc = m1 * (xc - mx1) + my1;
            }
            else
            {
                m1 = -(p2.X - p1.X) / (p2.Y - p1.Y);
                m2 = -(p3.X - p2.X) / (p3.Y - p2.Y);
                mx1 = (p1.X + p2.X) * 0.5;
                mx2 = (p2.X + p3.X) * 0.5;
                my1 = (p1.Y + p2.Y) * 0.5;
                my2 = (p2.Y + p3.Y) * 0.5;
                //Calculate CircumCircle center (xc,yc)
                xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
                yc = m1 * (xc - mx1) + my1;
            }

            double dx = p2.X - xc;
            double dy = p2.Y - yc;
            double rsqr = dx * dx + dy * dy;
            //double r = Math.Sqrt(rsqr); //Circumcircle radius
            dx = p.X - xc;
            dy = p.Y - yc;
            double drsqr = dx * dx + dy * dy;

            return (drsqr <= rsqr);
        }
    }
}
